/*
  File:    reinit.cpp
  Purpose: Utilities for reinitializing a level set function.  

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <iostream>

using namespace std;

#include "../ImprovedMath.h"
#include "../Matrix.h"
#include "WENO.h"
#include "extend.h"

double HeavisideDeltaX(double dXIn)
{
 extern double dx;
 extern double _PI;
 
 if( dXIn > dx )
 { return 1.0; }
 else if( dXIn < - dx )
 { return 0.0; }
 else
 { return (1.0/2.0)*(1.0 + dXIn / dx + (1.0/_PI)*sin(_PI*dXIn / dx)); }
}

// this one is the derivative of the numerical heaviside function

double HeavisideDeltaXDerivative(double dXIn)
{
 extern double dx;
 extern double _PI;
 
 if( fabs( dXIn) > dx )
 { return 0.0; }
 else
 { return (1.0/(2.0*dx))*(1.0+cos( dXIn * _PI / dx )); }
}

// this is the numerical sign function
double SignDeltaX(double dIn)
{ return 2.0 * ( HeavisideDeltaX(dIn) - 1.0/2.0); }

// uses third-order RK and 3rd or 5th-order WENO

double Output_Radius(Matrix &mIn)
{
 extern int Xnodes;
 extern int Ynodes;
 double output = 0.0;
 extern double _PI;
 
 for(int i=0;i<Xnodes;i++)
 {
	for(int j=0;j<Ynodes;j++)
	{
	
	if(*mIn(i,j)<0)
	{
	output++;
	}
	
	}	
 }
 output = output/_PI;
 output = sqrt(output)/100;//unit in mm
 
 return output;
}

Matrix Reinit(Matrix &mIn)
{
 int i=0, j=0;
 extern Matrix X;
 extern Matrix Y;
 extern double dx;
 extern double dy;
 extern int Xnodes;
 extern int Ynodes;
 extern int UpdateNodes;
 
 Matrix mOut = extend(mIn, 3);
 
 int M = mOut.TellRows();
 int N = mOut.TellCols();

 Matrix mSign = zeros(M,N);

 for(i=0; i < M ; i++)
 {
  for(j=0; j < N ; j++)
  { 
//   if( *mOut(i,j) < 3 ) // new on 11-19-2006 for bandedness
//   {
    *mSign(i,j) = SignDeltaX( *mOut(i,j) );
//   }
  }
 }
 
 Matrix mNext = mOut;

 Matrix mPrev;
 
 Matrix mOperator0;
 Matrix mOperator1;
 Matrix mOperator2;
 
// Matrix mNext0;
// Matrix mNext1;
 
 double dDeltaTtemp = .5 * dx;
 
 int k;
 for(k=1; k <= 2*UpdateNodes ; k++)
 {
  mPrev = mNext;
  mOperator0 = mSign - WENO5(mSign, mPrev); 
//  mOperator0 = mSign - WENO5banded(mSign, mPrev, 3); 
  mNext = mPrev + dDeltaTtemp * mOperator0;

  mOperator1 = mSign - WENO5(mSign, mNext);
//  mOperator1 = mSign - WENO5banded(mSign, mNext,3);
  mNext = mPrev + .25 * dDeltaTtemp * ( mOperator0 + mOperator1 );
  
  mOperator2 = mSign - WENO5(mSign, mNext);
//  mOperator2 = mSign - WENO5banded(mSign,mNext,3);
  mNext = mPrev +  (dDeltaTtemp/6.0)*( mOperator0 + mOperator1 + 4.0*mOperator2);
 }
 
 mOut = contract( mNext, 3);
 return mOut;
}
